home *** CD-ROM | disk | FTP | other *** search
/ Nebula 2 / Nebula Two.iso / SourceCode / MiscKit1.7.1 / MiscKit / Source / MiscGISKit / MiscPlanetCoordConverter.m < prev    next >
Text File  |  1995-07-20  |  11KB  |  351 lines
  1. /*====================== MiscPlanetCoordConverter.m =========================*/
  2. /* MiscPlanetCoordConverter class supports the calculations required for
  3.    converting between World and Universal Transverse Mercator coordinates.
  4.  
  5.    There is only one instance ever, so unless changes are made, this class
  6.    is NON REENTRANT.
  7.  
  8.    For information on the underlying mathematics, refer to:
  9.  
  10.      1- Ordinance Survey Information, "Transverse Mercator Projection,
  11.         Constants, Formula and Methods", March 1983
  12.  
  13.      2- Ordinance Survey, "Tables for the Transverse Mercator Projection
  14.         of Ireland", 1953, reprinted 1971
  15.  
  16. IMPORTANT: These equations are accurate to within 1 millimeter. Calculations
  17.        requiring greater accuracy must use formulae in:
  18.  
  19.         Redfern, JCB, "Transverse Mercator Formulae", 1948,
  20.         Empire Survey Review, 9(69) pg318-322
  21.  
  22.         Extra decimal places are stored only for the purpose of 
  23.         slowing error propagation that affects the numbers at the
  24.         millimeter scale, not because they are meaningful in and of 
  25.         themselves.
  26.  
  27.    DMA Release 0.8, Copyright @1993 by Genesis Project, Ltd. All Rights
  28.    Reserved. For further information on terms and conditions see:
  29.         Documentation/GISKit/Agreements-Legal-README
  30.  
  31. HISTORY
  32. 12-Mar-93  Dale Amon at GPL
  33.        Split UTMGrid into constants and converter parts.
  34. 22-Feb-93  Dale Amon at GPL
  35.        Created.
  36. */
  37.  
  38. #import <math.h>
  39. #import <misckit/miscgiskit.h>
  40.  
  41. @implementation MiscPlanetCoordConverter
  42.  
  43. /*===========================================================================*/
  44. /* Internal methods */
  45. /*===========================================================================*/
  46. /* calculations required for both directions of conversion. Caches values for
  47.    later use.
  48.  
  49.    Uses UTMConstants:
  50.     eSqrd = eccentricity squared
  51.     aF0   = major semiaxis
  52.  
  53. */
  54.  
  55. - (void) blackboardCalc: (double) phi
  56.   {    double    tmp;
  57.  
  58.     sinPhi  = sin(phi);
  59.     sin2Phi = sinPhi * sinPhi;
  60.  
  61.     tmp    = 1.0 - xlate->eSqrd * sin2Phi;
  62.     nu     =               xlate->aF0 / sqrt(tmp);
  63.     rho    = (nu * (1.0 - xlate->eSqrd)) / tmp;
  64.     etaSqrd= (nu/rho) - 1.0;
  65.  
  66.     return;
  67.   }
  68.  
  69.  
  70. /*---------------------------------------------------------------------------*/
  71. /* Calculate the Developed Arc of a Meridian from phi to the True Origin.
  72.    Uses local constants:
  73.     bF0  = semiminor axis, scaled.
  74.     phi0 = latitude of true origin
  75.    Caches values for later use.
  76. */
  77.  
  78. - (double) calcM: (double) phi
  79.   {    double    dif,sum;
  80.  
  81.     dif = phi - xlate->phi0;
  82.     sum = phi + xlate->phi0;
  83.  
  84.     M   = xlate->bF0 
  85.           * (xlate->M1 * dif - 
  86.              xlate->M2 * sin(    dif) * cos(    sum) +
  87.              xlate->M3 * sin(2.0*dif) * cos(2.0*sum) -
  88.              xlate->M4 * sin(3.0*dif) * cos(3.0*sum));
  89.  
  90.     return M;
  91.   }
  92.  
  93.  
  94. /*---------------------------------------------------------------------------*/
  95. /* Calculate our estimated latitide, phi prime.
  96.    Uses local constants:
  97.     aF0  = semimajor axis, acled
  98.     N0   = Grid north of True North
  99.     phi0 = latitude of true origin
  100.     convergence = difference at which we consider ourselves done
  101. */
  102.  
  103. - (double) calcPhiPrime: (double) N
  104.   {    double    tmp;
  105.     double    phiPrime = (N - xlate->N0)/xlate->aF0 + xlate->phi0;
  106.  
  107.     while(TRUE)
  108.      {    [self calcM: phiPrime];
  109.         tmp = N - xlate->N0 - M;
  110.         if (fabs(tmp) < xlate->convergence) break;
  111.         phiPrime += tmp/xlate->aF0;
  112.      }
  113.     return phiPrime;
  114.   }
  115.  
  116.  
  117. /*===========================================================================*/
  118. /* Conversion services */
  119. /*===========================================================================*/
  120. /* Calculate latitude and longitude given grid N and E. Accurate to 1 mm.
  121.    Uses constants:
  122.     lambda0 = longitude of grid origin
  123. */
  124.  
  125. - (void) utmToWorld
  126. {    double    E,N;
  127.     double    phiPrime;
  128.     double    y,y2,y3,y4,y5,y6,y7;
  129.     double    nu2,nu3,nu4,nu5,nu7;
  130.     double    tanPhiPrime,tan2PhiPrime,tan4PhiPrime,tan6PhiPrime;
  131.     double    secPhiPrime;
  132.     double    VII, VIII,IX, X, XI, XII, XIIA;
  133.  
  134.     /* conversion is driven by the source constants */
  135.     xlate = srcConstants;
  136.  
  137.     E = src[MISC_EASTINGS];
  138.     N = src[MISC_NORTHINGS];
  139.  
  140.     phiPrime = [self calcPhiPrime: N];
  141.  
  142.     /* precalculate nu, rho and etaSqrd and powers of nu */
  143.     [self blackboardCalc: phiPrime];
  144.     nu2   = nu*nu;
  145.     nu3   = nu2*nu;
  146.     nu4   = nu2*nu2;
  147.     nu5   = nu3*nu2;
  148.     nu7   = nu4*nu3;
  149.  
  150.     /* precalculate all the trig values */
  151.  
  152.     tanPhiPrime  = tan(phiPrime);
  153.     tan2PhiPrime = tanPhiPrime  * tanPhiPrime;
  154.     tan4PhiPrime = tan2PhiPrime * tan2PhiPrime;
  155.     tan6PhiPrime = tan4PhiPrime * tan2PhiPrime;
  156.     secPhiPrime  = 1/cos(phiPrime);
  157.  
  158.     /* precalculate all the powers of delta E */
  159.     y  = E - xlate->E0;
  160.     y2 = y*y;
  161.     y3 = y2*y;
  162.     y4 = y2*y2;
  163.     y5 = y3*y2;
  164.     y6 = y3*y3;
  165.     y7 = y4*y3;
  166.  
  167.     /* Calculate the terms used by the latitude and longitude equations */
  168.  
  169.     VII  = tanPhiPrime/(2.0*rho*nu);
  170.     VIII = tanPhiPrime/(24.0*rho*nu3) *
  171.         (5.0 + 3.0*tan2PhiPrime + etaSqrd - 9.0*etaSqrd*tan2PhiPrime);
  172.     IX   = tanPhiPrime/(720.0*rho*nu5) *
  173.         (61.0 + 90.0*tan2PhiPrime + 45.0*tan4PhiPrime);
  174.  
  175.     X    = secPhiPrime/nu;
  176.     XI   = secPhiPrime/(6.0*nu3) * (nu/rho + 2.0*tan2PhiPrime);
  177.     XII  = secPhiPrime/(120.0*nu5) *
  178.         (5.0 + 28.0*tan2PhiPrime + 24.0*tan4PhiPrime);
  179.     XIIA = secPhiPrime/(5040.0*nu7) *
  180.         (61.0 + 662.0*tan2PhiPrime +
  181.          1320.0*tan4PhiPrime + 720.0*tan6PhiPrime);
  182.  
  183.     /* and finally, we give you the latitude and the longitude */    
  184.     dst[MISC_LATITUDE]  = phiPrime     - y2*VII + y4*VIII + y6*IX;
  185.     dst[MISC_LONGITUDE] = xlate->lambda0 + y*X - y3*XI + y5*XII - y7*XIIA;
  186.     dst[MISC_ALTITUDE]  = src[MISC_ELEVATION];
  187. }
  188.  
  189.  
  190. /*---------------------------------------------------------------------------*/
  191. /* Given World Coordinates, calculate the local grid coordinates in
  192.    a UTM projection.
  193.    Uses constants:
  194.     N0,E0 = True origin offset from grid False origin
  195.     lambda0 = longitude of True origin. 
  196. */
  197.  
  198. - (void) worldToUTM
  199. {    double    phi,lambda;
  200.     double    cosPhi,cos2Phi,cos3Phi,cos5Phi;
  201.     double    tanPhi,tan2Phi,tan4Phi;
  202.     double    P1, P2, P3, P4, P5, P6;
  203.     double    I,II,III,IIIA,IV,V,VI;
  204.  
  205.     /* conversion is driven by the destination constants */
  206.     xlate = dstConstants;
  207.  
  208.     /* assumes the internal storage format is double precision radians */
  209.     phi = src[MISC_LATITUDE]; lambda = src[MISC_LONGITUDE];
  210.  
  211.     /* Calculate the Developed Arc of a Meridian from phi to the
  212.        True Origin. */
  213.     [self calcM: phi];
  214.  
  215.     /* calculate nu, rho and etaSqrd */
  216.     [self blackboardCalc: phi];
  217.  
  218.     /* precalculate all the trig values that are not on blackboard */
  219.     cosPhi  = cos(phi);
  220.     cos2Phi = cosPhi*cosPhi;
  221.     cos3Phi = cos2Phi*cosPhi;
  222.     cos5Phi = cos2Phi*cos3Phi;
  223.  
  224.     tanPhi  = tan(phi);
  225.     tan2Phi = tanPhi  * tanPhi;
  226.     tan4Phi = tan2Phi * tan2Phi;
  227.  
  228.     /* Precalculate all the powers of delta lambda */
  229.     P1 = lambda - xlate->lambda0;
  230.     P2 = P1 * P1;
  231.     P3 = P2 * P1;
  232.     P4 = P2 * P2;
  233.     P5 = P3 * P2;
  234.     P6 = P3 * P3;
  235.  
  236.     /* Calculate the terms used by the Easting and Northing equations */
  237.  
  238.     I   = M + xlate->N0;
  239.     II  = nu/2.0   * sinPhi * cosPhi; 
  240.     III = nu/24.0  * sinPhi * cos3Phi * (5.0 - tan2Phi + 9.0 * etaSqrd);
  241.     IIIA= nu/720.0 * sinPhi * cos5Phi * (61.0 - 58.0 * tan2Phi + tan4Phi);
  242.     IV  = nu       * cosPhi;
  243.     V   = nu/6.0   * cos3Phi * (nu/rho - tan2Phi);
  244.  
  245.     VI  = nu/120.0 * cos5Phi * 
  246.         (5.0 - 18.0*tan2Phi + tan4Phi + 14.0*etaSqrd - 
  247.          58.0 * etaSqrd * tan2Phi);
  248.  
  249.     /* And finally, we calculate the local grid coordinates */
  250.     dst[MISC_NORTHINGS] = I         + P2*II + P4*III + P6*IIIA;
  251.     dst[MISC_EASTINGS]  = xlate->E0 + P1*IV + P3*V   + P5*VI;
  252.     dst[MISC_ELEVATION] = src[MISC_ALTITUDE];
  253.  }
  254.  
  255.  
  256. /*===========================================================================*/
  257. /* Class methods */
  258. /*===========================================================================*/
  259. /* Initialize the class */
  260.  
  261. + initialize
  262.   {[MiscPlanetCoordConverter setVersion:MISC_PLANET_CONVERT_VERSION_ID]; return self;}
  263.  
  264. /*---------------------------------------------------------------------------*/
  265. /* Only one converter of this type is every needed. Of course if we got
  266.    into a really big multiprocess agora system there might be a case
  267.    for multiple convertors. Since this object is shared by many, it
  268.    doesn't particularly matter what zone it is allocated from.
  269. */
  270.  
  271. static id thePlanetCoordConverter = nil;
  272.  
  273. + new
  274. {
  275.     if (!thePlanetCoordConverter)
  276.     thePlanetCoordConverter = [[super alloc] init];
  277.     return thePlanetCoordConverter;
  278. }
  279.  
  280.  
  281. /*---------------------------------------------------------------------------*/
  282. /* Since there is only one object needed, alloc and allocFromZone: are 
  283.    disabled and considered errors. free is simply overridden and made into a
  284.    no op.
  285.  
  286.    Superalloc is included because we want to allow our subclasses to also 
  287.    create single private instances.
  288. */
  289.  
  290. +alloc 
  291.  {    [self error:"  %s class should not be sent '%s' messages\n",
  292.             [[self class] name], sel_getName(_cmd)];
  293.     return self;
  294.  }
  295.  
  296. +allocFromZone: (NXZone *)zone
  297.  {    [self error:"  %s class should not be sent '%s' messages\n",
  298.             [[self class] name], sel_getName(_cmd)];
  299.     return self;
  300.  }
  301.  
  302. - free {return self;}
  303.  
  304. +superalloc {return [super alloc];}
  305.  
  306. /*===========================================================================*/
  307. /* Initialization methods */
  308. /*===========================================================================*/
  309. /* We add a list of all the services which we are able to provide.
  310.    Note that there is only one PlanetCoordConverter ever created. Once 
  311.    initialized the same object is given to all comers and can not be destroyed.
  312. */
  313.  
  314. - init
  315.   { Class    world, utm, ukutm, zutm;
  316.  
  317.     [super init];
  318.     world  = [MiscWorldCoord   class];
  319.     utm    = [MiscUTMCoord     class];
  320.     ukutm  = [MiscUKUTMCoord   class];
  321.     zutm   = [MiscZoneUTMCoord class];
  322.  
  323.     [self addService: @selector(utmToWorld)   convertsFrom: zutm  to: world]; 
  324.     [self addService: @selector(worldToUTM)   convertsFrom: world to:  zutm]; 
  325.     [self addService: @selector(utmToWorld)   convertsFrom: ukutm to: world]; 
  326.     [self addService: @selector(worldToUTM)   convertsFrom: world to: ukutm]; 
  327.     [self addService: @selector(utmToWorld)   convertsFrom: utm   to: world]; 
  328.     [self addService: @selector(worldToUTM)   convertsFrom: world to: utm]; 
  329.     [self addService: [self fastCopySelector] convertsFrom: ukutm to: ukutm]; 
  330.     [self addService: [self fastCopySelector] convertsFrom: world to: world]; 
  331.  
  332.     /* This will method return NO if they are in different zones */
  333.     [self addService: [self fastCopySelector] convertsFrom: zutm to: zutm]; 
  334.  
  335.     return self;
  336.   }
  337.  
  338.  
  339. /*===========================================================================*/
  340. /* Archive methods */
  341. /*===========================================================================*/
  342.  
  343. - finishUnarchiving
  344. {
  345.     [self free];
  346.     return [MiscPlanetCoordConverter new];
  347. }
  348.  
  349.  
  350. @end
  351.